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stability  of  the  guidance/control  scheme  is  analyzed.  Results  are  presented  based  on  time 
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I .   INTRODUCTION 

A .   BACKGROUND 

Marine  vehicles  that  require  high  maneuverability,  such  as 
hydrofoils,  require  quick  response,  control  and  guidance  to 
maintain  accurate  track  keeping  characteristics.  This  is 
accomplished  through  successful  path  planning,  navigation, 
guidance  and  autopilot  design  [Ref .  1] . 

Sufficient  information  is  obtained  from  charted  obstacles 
and  the  environment  for  smooth  paths  to  be  generated  for  the 
vehicle  to  follow  [Ref.  2]  .  A  certain  level  of  feedback  is 
provided  through  the  use  of  sonar  and  acoustics  in  order  to 
replan  a  path  when  uncharted  obstacles  are  encountered  or  when 
the  mission  objectives  are  changed.  Based  on  inertial 
positional  information,  the  guidance  law  provides  the 
appropriate  vehicle  heading  by  suitable  use  of  the  vehicle 
effectors  such  as  rudders,  dive  planes,  and  cross  body 
thrusters.  However,  lags  in  obtaining  and  processing  inertial 
positional  information  can  limit  vehicle  reaction  time.  In 
addition,  the  guidance  law  must  be  as  fast  as  possible  in 
order  to  ensure  accurate  path  keeping  characteristics  [Ref. 
3]  .  Therefore,  the  stability  of  the  combined  guidance/control 
scheme  becomes  an  issue  that  needs  to  be  analyzed. 


B.  OBJECTIVE  OF  THIS  THESIS 

This  thesis  investigates  the  effects  of  positional 
information  time  lags  on  guidance  and  control.  The 
relationship  of  the  critical  visibility  (minimum  vehicle 
lookahead  distance)  versus  the  controller  time  constant  and 
its  effect  on  the  stability  of  the  combined  guidance/control 
scheme  is  analyzed.  Results  are  presented  based  on 
computations  using  time  domain  and  frequency  response 
techniques.  All  computations  are  performed  for  a  dynamic 
model  of  the  NPS  AUV  II  [Refs.  4  and  5]  for  which  a  complete 
set  of  hydrodynamic  coefficients  and  geometric  properties  is 
available . 

C.  THESIS  OUTLINE 

Chapter  II  presents  the  vehicle  equations  of  motion  in  the 
horizontal  plane,  the  equations  that  govern  steering  control, 
and  the  guidance  scheme  used  in  this  research. 

Chapter  III  presents  the  stability  analysis  based  on 
eigenvalue  computation  of  the  combined  guidance  and  control 
scheme  developed  in  Chapter  II.  Results  are  presented  based 
on  the  stability  curves  developed  for  first,  second  and  third 
order  approximations  for  time  lag  in  the  commanded  vehicle 
heading  in  the  control  law. 

Chapter  IV  presents  an  exact  computation  of  the  stability 
curves  based  on  frequency  response  methods,  which  utilizes  the 
Nyquist  criterion  for  stability. 


II.  EQUATIONS  OF  MOTION,  CONTROLLER  DESIGN  AND  GUIDANCE  SCHEME 

A.  INTRODUCTION 

The  equations  of  motion  that  describe  maneuvering  of  a 
marine  vehicle  in  the  horizontal  plane  are  presented  in  this 
chapter.  A  linear  full  state  steering  feedback  control  law  is 
designed  based  on  the  three  equations  in  sway,  yaw,  and  rate 
of  change  of  heading  angle.  The  guidance  scheme  is  based  on 
pure  pursuit  guidance  for  path  following  along  straight  line 
segments . 

B.  EQUATIONS  OF  MOTION 

Restricting  our  attention  to  motion  in  the  horizontal 
plane  (steering  control),  the  mathematical  model  consists  of 
the  nonlinear  sway  and  yaw  equations  of  motion  only.  With  a 
moving  coordinate  frame  fixed  at  the  vehicle's  geometric 
center,  Newton's  equations  of  motion  are 

m( v+ur+xGf-yGr2)  -Y  (2.1) 

Izr  +mxG  {v+ur)  -myGvr =N  (2.2) 


where  v  and  r  are  the  relative  sway  and  yaw  velocities  of  the 
moving  vehicle  with  respect  to  the  water;  m  is  the  mass  of  the 
vehicle;  Xq,  yG  are  the  respective  lateral  and  longitudinal 
positions  of  the  center  of  gravity;  and  Y,N  represent  the 
total  excitation  sway  force  and  yaw  moment,  respectively. 
These  forces  can  be  expressed  as  the  sum  of  quadratic  drag 
terms  and  first  order  memoryless  polynomials  in  v  and  r,  which 
represent  the  added  mass  and  damping  due  to  the  vehicle's 
motion  through  the  water.  In  this  way,  Y  and  N  can  be 
represented  by 


Y  -  ■P-PY,f  +  £-P(Y^v+Yrur)+£-t2Yvuv 
2^2     v         r  2        v 


2  J  tan    °y  v+Er      2   6 


N  =  -P  fi^  r  +  _P  <>4  (N  v+Nriir)  +  -P-PNvuv 
2        *    2     v         r  2        v 

2  J  tan    "y  \v+lr\  2      ° 

where  C  is  the  vehicle  length,  and  Ya,  Na  represent  partial 
derivatives  of  Y  and  N  with  respect  to  a,  C^  is  the  drag 
coefficient,  and  8  is  the  rudder  angle. 

Equations  (2.1)  and  (2.2)  can  be  nondimensionalized  with 
respect  to  the  constant  forward  speed  u,  and  the  vehicle 
length  i  ;  with  the  dimensionless  time  variable  being  tu/d  . 


The  cross  flow  integral  drag  terms  become  important  for 
hovering  operations  or  low  speed  manuevering,  whereas  at 
relatively  high  speeds  and  low  angles  of  attack  (with  respect 
to  the  water),  the  steering  response  is  predominantly  linear. 
The  model  becomes  complete  with  the  addition  of  the 
expressions  for  the  vehicle  yaw  rate 

i|r=r  (2.3) 

and  the  inertial  position  rate  (horizontal  plane) 

y=usini|r  +  vcosi|r  (2.4) 

where  \|/  is  the  heading  angle  as  shown  in  Figure  1. 

Equations  (2.1),  (2.2),  and  (2.3)  can  be  written  as  a  set 
of  three  coupled  linear  differential  equations  of  the  form 

t|/=r  (2.5a) 


v^=a11uv+a12ur+jb1u26  (2.5b! 


r=al2uv+a2Zur+b2u26  (2.5c 
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Figure  1.   Vehicle  Geometry  and  Definition  of  Symbols 


where  the  coefficients  ai;j  and  b1  are  functions  of  the  vehicle 
hydrodynamic  and  geometric  properties,  \\f  is  the  heading  angle 
and  5  is  the  rudder  angle. 

The  linearized  form  of  equation  (2.4)  at  a  nominal  u  is 

y=ui|;  +  v  (2.6) 

Hence,  the  final  set  of  linearized  equations  of  motion  used 
for  steering  control  are 

i|/=r  (2.7a) 

^=a11uvr+a12ur+jb1u25  (2.7b) 


f=a12uv+a22ur+b2u26  (2.7c) 


y=wty  +  v  (2 .7d) 

at  any  nominal  u. 

This  system  of  equations  becomes  the  basis  for  the 
controller  design. 


C.   CONTROLLER  DESIGN 

Equations  (2.7a),  (2.7b),  and  (2.7c)  govern  the  steering 
control  of  the  model  used  in  this  research. 

These  equations  can  be  represented  in  the  form 

x=Ax+bb  (2.8) 

where  the  state  vector  equation  is 

x   =  [i|r,  v,r]T  (2.9) 

Linear  full  state  feedback  is  introduced  in  the  form 

fi=/c1(T|f-i|fc)  +k2v+k3r  (2.10) 

where  \|/c  is  the  commanded  heading  angle,  and  the  gains  kw  k2, 
k3  are  computed  by  pole  placement  such  that  the  closed  loop 
system  of  equations  (2.8),  (2.10)  has  the  desired  dynamics. 
The  existence  of  the  difference  (\|/-\j/c)  in  the  control  law 
(2.9)  has  the  effect  of  trying  to  point  the  vehicle's 
longitudinal  axis  towards  the  desired  heading. 

If  the  desired  characteristic  equation  is  selected  so  that 
its  time  constant  is  tc  dimensionless  seconds,  its  general 
form  becomes 


X3+axA.2+a2A.+a3  =0  (2.11 


where 


3  3  1 

C  U  c  L  c 


The  controller  gains  are  computed  from 


kx= ^ (2.12a 

{b2axl-bra21)  u3 


(jb^^ -Jb2a12)  u3^2+(jb2a11-jb1a21)u3ic3=a2+iD2u2ic1       (2.12b) 

Jb1u2ic2+jb2u2ic3=-a1-(a11+a22)ii  (2.12c) 

Due  to  the  explicit  dependence  on  u,  the  gains  k1(  k2,  k3  are 
continuously  updated  for  different  nominal  forward  speeds. 

D.   GUIDANCE  SCHEME 

The   guidance   scheme   used   in   this   research   is   an 
orientation  based  command  scheme.    In  this   scheme,   the 


commanded  vehicle  heading  angle  \|/-  is  a  function  of  the  line 
of  sight  angle  a  between  the  actual  vehicle  position  and  a 
reference  position  on  the  nominal  path  which  is  located  at  a 
constant  lookahead  or  visibility  distance  d  ahead  of  the 
vehicle  as  shown  in  Figure  1. 

This  line  of  sight  angle  is  defined  by 


tano  =  -^  (2.13) 

d 


The  autopilot  is  then  called  upon  to  deliver  the  commanded 
heading  angle  \|/c .   The  simplest  orientation  guidance  law  is 
pure  pursuit  guidance  where  the  commanded  heading  angle  \|/c  in 
the  control  law  (2.9)  equals  the  line  of  sight  angle  a  in 
(2.13) . 


Then  \j/r  is  defined  by 


ili^  =  -arctan^,  (2  .14 

c  d 


For  relatively  small  angles 


i|r_  =  -arctan^  «  -i£  (2.15) 


10 


This  results  in  a  three  state  heading  autopilot  design 
which  can  exhibit  very  robust  characteristics.  However,  the 
stability  of  this  scheme  may  become  an  issue  when  transiting 
along  straight  line  segments  since  the  commanded  heading  angle 
is  a  function  of  the  vehicle  response,  i.e.,  autopilot 
response  is  limited  by  time  lags  in  vehicle  dynamics.  We  will 
next  examine  the  stability  of  this  scheme. 


11 


III.   EIGENVALUE  ANALYSIS 

A.  INTRODUCTION 

In  this  chapter,  we  present  the  stability  analysis,  based 
on  eigenvalue  computation,  of  the  combined  guidance  and 
control  scheme  developed  in  Chapter  II.  A  brief  presentation 
of  the  stability  curve  is  developed  first  for  the  case  of  zero 
time  lag.  Incorporation  of  a  nonzero  time  lag  in  the 
positional  information  formally  results  in  an  infinite  state 
space  system.  This  is  truncated  into  first,  second,  and  third 
order  approximations.  The  resulting  eigenvalue  problems  are 
then  analyzed  with  emphasis  on  stability  curves. 

B.  STABILITY  CONSIDERATIONS  AND  TIME  LAG 

1.   Stability  Considerations 

In  pure  pursuit  guidance,  where  the  commanded  heading 
angle  \j/c  in  the  control  law  equals  the  line  of  sight  angle  a, 
the  trivial  equilibrium  state  which  corresponds  to  straight 
line  motion  is  characterized  by 

ij;  =  v=r=y=0  (3.1) 

Linearization  of  the  state  equations  in  the  vicinity  of  (3.1) 
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produces  the  linear  system 


x=Ax,  (3.2) 

where  the  complete  state  vector  is 

x=[i|r,  v,r,y]T  (3.3) 

Local  stability  properties  of  (3.1)  are  then  established  by 
the  eigenvalues  of  A.  The  characteristic  equation  of  (3.2)  is 
a  quartic  of  the  form 

\*+B\2+Ck2+Dk+E=Q,  (3.4) 

where  the  coefficients  B,  C,  D,  E  are  functions  of  the  vehicle 
properties,  the  lookahead  distance  d,  and  the  controller  gains 
klt    k2,  k3 .   A  pair  of  complex  conjugate  roots  of  (3.4)  crosses 
the  imaginary  axis  when  the  following  conditions  are  met 
[Ref .  6] . 

BCD-B2E-D2=Q,  (3.5a) 


->0  (3.5b! 

B 
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These    two    conditions    (3.5a)    and    (3.5b)    translate    to 


a1d2+a2d+a3=0,  (3.6a) 


d>hra22-b2a12-b  (3>6b) 


■^2ail     -^Ia21 


where 


a1=g1g2-a3  (3 .7a) 


a  _  (g1a2-2a3)  (iy322-.b2a12-.b2)  _         ^gigj         _g2u 
2  ^auAa2i  (^^ii -^21) u      1     (3-7b: 


a  _  -  (^ia22 --b2a12--b2)  [,b1tt1+(,b1a22-,b2a12-,b2)u]tt3  ?c 

(jhaa11-Jb1a21)2u 


The  conditions  (3.5a)  and  (3.5b)  determine  the  critical 
visibility  dcrit  for  stability.  For  d  >  dcrit/  the  equilibrium 
state  (3.1)  is  stable  which  means  that  the  control  law  will 
drive  and  keep  the  vehicle  onto  the  straight  line  path.  For 
d  <  dcric/  the  equilibrium  state  loses  its  stability  and  the 
vehicle  response  becomes  oscillatory  as  a  result  of  the  pair 
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of  complex  conjugate  eigenvalues  with  positive  real  parts. 
Results  of  the  critical  visibility  versus  the 
controller  time  constant  tc  for  zero  time  lag  are  presented  in 
Figure  2.  As  expected,  the  higher  values  of  tc  (i.e., softer 
controller)  require  higher  lookahead  distances  d  for  path 
accuracy.  This  agrees  with  existing  guidelines  that  the 
guidance  law  must  be  sufficiently  slower  than  the  controller 
in  order  for  the  dynamics  of  the  two  not  to  interfere  with 
each  other  and,  therefore,  guarantee  stability  of  the  combined 
scheme.  Very  high  values  of  d  correspond  to  a  very  slow 
guidance  law  with  a  loss  in  speed  of  response  and  path 
accuracy  [Ref .  1] . 

2 .   Time  Lag 

In  ocean  vehicles,  all  states  needed  in  the  control 
law  are  readily  available  at  the  required  rate,  with  the 
possible  exception  of  the  positional  information  y  (Figure  1) . 
The  latter  may  require  a  significant  data  analysis  and 
reduction  of  sonar  returns  and  inertial  navigational 
information.  As  a  result,  it  is  likely  that  a  time  lag  will 
exist  in  the  positional  information  y  which  is  used  in  the 
guidance  law. 

With  the  introduction  of  a  time  lag  T  (sec)  ,  the  commanded 
rudder  angle  \|/c  in  the  pursuit  guidance  control  law  becomes 
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Figure  2.   Critical  Visibility  Distance  d  versus  tc  for 
Zero  Time  Lag  T. 
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♦„  =  -arctanil^ 
c  d 


and  for  relatively  small  angles 


t  ..zJi-n  (3.9) 


The  resultant  rudder  control  law  (2.9)  becomes 


5=ic1\|f+Jc1y(^ ,  T)  +k2v+k3r  (3.10 


After  some  algebra,  the  resultant  linearized  equations  of 
motion  (2.7a),  (2.7b),  (2.7c),  (2.7d)  become 


i|r=r  (3.11a 


v=b1u2ic1!|r  +  (a11u+i31u2^:2)  v+  {a12u+bxuzk2)  r+bxu2-£y(  t-T)     (3  .lib; 


i-=jb2u2ic1i|r+(a21Lz+i32Lz2/c2)  v+  (a22u+-b2u2£3)  r+Jb2u2-^y  ( t-T)     (3  .lie 


17 


y=uty  +  v  (3. lid! 


The  Taylor  expansion  for  the  term  y(t-T)  is 


y(t-T)  *  y-Ty+  —  y-  —  y+.  .  .  (3.12) 

2    6 


We  can  now  examine  the  stability  of  the  previously  developed 
control  scheme  for  first,  second  and  third  order 
approximations  for  time  lag  in  the  control  law. 

C.   FIRST  ORDER  APPROXIMATION  FOR  TIME  LAG 

For  a  first  order  approximation  for  time  lag  T 


y(t-T)  =  y-Ty  (3  .13) 


where 


y=m|/  +  v  (3.14) 

After  some  algebra,  the  resultant  linearized  equations  of 
motion  (3.11a),  (3.11b),  (3.11c),  (3. lid)  become 


i|r=r 


3 .15a) 


*i 


K 


^=(i?1u2ic1-Jb1u3-^D*+(a11u+i)1u2^2-2)1u2-^D  v+(a12u+jb1u2/c3)r+2)1u2-^y(3  >15b 


K 


»*i 


2   ^1. 


r=  (b^k^u^T)  ty+  (a21u+b2u2k2-b2u2^T)  v+  (a22u+b2u2k2)  r+b2u2-+y{3    1Sc 


y=ui|i  +  v 


3.15di 


The  state  space  form  of  the  equations  of  motion  is 


x=Ax 


(3  .16 


In  matrix  form,  (3.15a),  (3.15b),  (3.15c),  (3.15d)  become 


(JhluaJcl-2xlu,iT)    (a^u+bxuzk2-b^u2  ^T)    {a^u+b^k^)     (2)tu2-^) 
(ij2u2^-fc2u^D    (a21u+i2u2Jc2-i32u2^D    (a^u^u'*,)     (£>2"2:§> 


(3.17) 


The  standard  eigenvalue  problem  (3.17)  is  solved  to 
determine  numerically  the  critical  d  versus  tc  curve  for  a 
given  value  of  T.   Appendix  A  presents  the  program  used  for 


19 


the  first  order  approximation  for  time  lag. 

Figure  3  shows  the  resulting  stability  curves  for  the 
first  order  case.  It  shows  the  relation  between  the  critical 
visibility  d  and  the  controller  time  constant  tc  for  values  of 
time  lag  of  0.0,  0.5,  1.0,  1.5  and  2.0  seconds. 

As  expected,  the  stability  curve  shifts  to  the  right  with 
increasing  values  of  time  lag.   For  the  same  time  constant  tc, 
the  critical  visibility  distance  increases  with  corresponding 
increases  in  the  amount  of  time  lag. 

D.   SECOND  ORDER  APPROXIMATION  FOR  TIME  LAG 

For  a  second  order  approximation  for  time  lag  T,  we  have 


T2 


y(  t-T)  =  y-Ty+—y  (3.1: 


where 


y=uty  +  v,     y^ur+v  (3.14),  (3.19) 

Following  some  algebra,  the  resultant  linearized  equations  of 
motion  (3.11a),  (3.11b),  (3.11c),  (3. lid)  become 
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Figure  3.   First  Order  Approximation  For  Time  Lag: 

Critical  Visibility  Distance  d  versus  tc  for 
Time  Lags  T  =  0.0,  0.5,  1.0,  1.5,  2.0  sec. 
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where 


i|r  =  r 


(3  .20a 


^=A2 1  ^ +^2  2  ^+^2  3  r  +  A>4y 


3.20b) 


f=A3iq+A32v+A33r+A34y 


(3 .20c) 


y-u^+v 


(3.20d) 


*21 


1      2d 


x22 


a^u^u2^-^2-^ 

1  -2),  u 2  — ^  r2 

1      2d 


a ,  o  u  +b,  u  2 ic,  +  b,  u 3  — \  Tz 


'12"    ~1"    ix3     ~1' 


A 


2d 
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i-£nu2  — r2 
1      2d 


22 


b,u2^ 
A 


24  k 

l-b^—^T2 


2d 


,3*1 


k  k  (b^k.-b^-^T) 

A31  =    (b2u*k,-b2u^T)+{b2u*-l-T2) 


^  ^^ 


.  .  {a^u+b^uzk2-biU2-±T) 

A22  =   (a21u+b2u2k2-b2u2-±T^(b2u*-j-T2) 


k 
( a12  u  +jbx  u  "k2  +jbx  u J 

A,,  =   (a„u+2?,u2Jt,+ij,u3-^i7'2)  +  (£,u2— ,T2) 


2d  2      2d  Jq 


(a12u+2>1ualc3+Jb1u3-^:r2) 


2d 


a34  =  (h.u^J  +  ^u'-^r') 


(l-i^u2—^2) 
1      2d 


As  with  the  first  order  case,  the  state  space  equations  of 
motion  reduce  to 
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k-Ax 


3.16 


In  matrix  form,  (3.20a),  (3.20b),  (3.20c),  (3.20d)  become 


* 

V 

- 

r 

y 

■"21    "^22    "^23    ^24 


An    ^32    ^33     A 


34 


* 


V 


y 


(3.21) 


Simular  to  the  first  order  case,  the  standard  eigenvalue 
problem  (3.21)  is  solved  using  Fortran  programming  and  MATLAB 
in  order  to  develop  the  stability  characteristics  of  the 
second  order  case.  Appendix  B  presents  the  program  used  to 
develop  the  stability  characteristics  for  this  case. 

Figure  4  shows  the  resulting  stability  curves  for  this 
case.   The  results  of  the  critical  visibility  d  versus  the 
controller  time  constant  tc  for  time  lags  of  0.0,  0.5,  1.0, 
1.5  and  2.0  sees  are  shown. 

As  with  the  first  order  case,  the  higher  values  of  tc 
require  higher  lookahead  distances  d  for  path  accuracy.  High 
values  of  d  correspond  to  a  slow  guidance  law  with  a  loss  in 
speed  of  response  and  path  accuracy,  and  the  stability  curve 
shifts  to  the  right  with  increasing  values  of  time  lag. 
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Figure  4.   Second  Order  Approximation  For  Time  Lag: 

Critical  Visibility  Distance  d  versus  tc  for 
Time  Lags  T  =  0.0,  0.5,  1.0,  1.5,  2.0  sec. 


25 


However,  due  to  being  more  accurate  than  the  first  order  case, 
the  second  order  approximation  provides  tighter  stability 
curves  for  the  same  time  lags. 

E.   THIRD  ORDER  APPROXIMATION  FOR  TIME  LAG 

For  a  third  order  approximation  for  time  lag  T 


y(t-T)  =  y-Ty+^y-Tly  (3.22 

2     6 


where 


y=uty+v,     y=ur+v,     y-uz^v        (3.14),  (3.19),  (3.23) 

After  a  considerable  amount  of  algebra  and  following  a  simular 
procedure  as  with  both  the  previous  two  cases,  the  resultant 
linearized  equations  of  motion  reduce  to  the  state  space  form 

Bx=Ax  (3.24) 


or 


x=B~1Ax 
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The  resultant  linearized  equations  of  motion  (3.11a), 
3.11b),  (3.11c),  (3. lid)  become 

ijr=r  (3.25a) 

B33±+B3sV2  =  A^y+A^v^A^Vz+A^r+A^y  (3.25b) 

B53^+fl55*2     =   A5iy+A52Vl+A52V2+A5ir+ASSy  (3.25c) 

y=ui^+v1  (3.25d) 


where 


vx  =  v2  ( 3  .  2  5  e ) 


33         x       6d 


s35  =  W  4r3 

35         1       6d 


A31  =  b^k^b^-^T 
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and 


A>2  =  axxu^bxu2k2-bxu2-^T 


^33  =  ax2u+bxu2k3+bxu^-^T2 


^34  ^1U     ^j 


35        x      2d 


53  2       6d 


S55  =  *2"2^3 


A5l  =  b2u2-b2u>-±T 


Asz  =  a21u+i?2u2ic2-2)2u2-^r 
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^53    =    Z22U+t>2U2k2+b2U2^T2 


^54  =  £2"2^ 


A„  =  b?u2^-T2 
55    2   2d 


In  matrix  form,  (3.25a),  (3.25b),  (3.25c),  (3.25d),  (3.25e) 
become 


0        10         0         0 
0       0      B33       0      s35 
0        0         0         10 


0        0       B53       0        B 


55 


* 


0  0  0  0 


-^31  "^32  "^33  Au        Ai5 


■^51  -^52  ^53  "^54         ^*55 


* 
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Simular  to  both  the  first  order  and  second  order  cases, 
the  generalized  eigenvalue  problem  (3.21)  is  solved  using 
Fortran  programming  and  MATLAB  in  order  to  develop  the 
stability  characteristics  for  the  third  order  case.  Appendix 
C  presents  the  program  used  to  develop  the  stability 
characteristics  for  this  case. 

Figure  5  shows  the  resulting  stability  curves  for  this 
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Figure  5.   Third  Order  Approximation  For  Time  Lag: 

Critical  Visibility  Distance  d  versus  tc  for 
Time  Lags  T  =  0.0,  0.5,  1.0,  1.5,  2.0  sec. 
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case.  The  results  of  the  critical  visibility  d  versus  the 
controller  time  constant  tc  for  time  lags  of  0.0,  0.5,  1.0, 
1.5  and  2.0  seconds  are  shown.  As  with  the  first  and  second 
order  cases,  the  higher  values  of  tc  require  higher  lookahead 
distances  d  for  path  accuracy,  high  values  of  d  correspond  to 
a  slow  guidance  law  with  a  loss  in  speed  of  response  and  path 
accuracy,  and  the  stability  curve  shifts  to  the  right  with 
increasing  values  of  time  lag.  However,  •  due  to  being  more 
accurate  than  either  of  the  other  two  cases,  the  third  order 
approximation  provides  the  tightest  stability  curves  for  the 
same  time  lags. 

F.   COMPARISON  OF  RESULTS 

Results  of  the  critical  visibility  d  versus  the  controller 
time  constant  tc  for  first,  second  and  third  order 
approximations  for  time  lag  have  been  obtained.  Figures  3,  4 
and  5  presented  the  stability  curves  for  first,  second  and 
third  order  approximations  for  time  lags  of  0.0,  0.5,  1.0,  1.5 
and  2.0  sees.  These  figures  have  shown  that  stability  of  the 
control  scheme  decreases  with  increasing  time  lag,  and  for  the 
same  time  constant,  the  critical  visibility  distance  increases 
with  increasing  time  lag. 

Figures  6,  7,  8  and  9  present  a  comparison  of  first, 
second  and  third  order  approximations  for  each  time  lag 
considered.  These  figures  show  that  for  each  time  lag,  the 
third  order  model  presents  the  largest  region  of  stability  for 
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Figure  6.   Critical  Visibility  Distance  d  versus  tc  for 
First,  Second  and  Third  Order  Approximation 
Cases  and  Time  Lag  T  =  0 . 5  sec . 
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Figure  7.   Critical  Visibility  Distance  d  versus  tc  for 
First,  Second  and  Third  Order  Approximation 
Cases  and  Time  Lag  T  =  1 . 0  sec . 
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Figure  8.   Critical  Visibility  Distance  d  versus  tc  for 
First,  Second  and  Third  Order  Approximation 
Cases  and  Time  Lag  T  =  1.5  sec. 
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Figure  9.   Critical  Visibility  Distance  d  versus  tc  for 
First,  Second  and  Third  Order  Approximation 
Cases  and  Time  Lag  T  =  2 . 0  sec . 
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vehicle  straight  line  path  accuracy.  In  addition,  for  the 
same  time  constant,  the  critical  visibility  is  least  for  the 
third  order  case.  This  demonstrates  that  the  minimum 
lookahead  distance  required  for  straight  line  motion  stability 
decreases  with  controller  time  lag  accuracy. 

It  can  be  seen  that  the  first  order  approximation  is  the 
most  conservative  for  design  purposes  since  it  predicts  the 
highest  required  minimum  visibility  distance.  As  expected, 
the  differences  among  the  three  approximations  are  more 
pronounced  as  the  time  lag  increases  and  as  the  controller 
time  constant  becomes  smaller;  i.e.,  tighter  control. 

Figures  10,  11,  12  and  13  present  the  differences  in 
critical  visibilities  among  first,  second  and  third  order 
approximations  for  each  individual  time  lag  considered.  A 
comparison  of  these  figures  demonstrates  that  the  difference 
in  critical  visibility  between  respective  orders  (time  lag 
models)  decreases  with  decreasing  time  lag.  This  indicates 
that  for  low  time  lags,  controller  time  lag  accuracy  is  not  as 
crucial  for  maintaining  vehicle  straight  line  path  accuracy  as 
that  of  higher  time  lags.  The  greater  the  time  lag,  the 
greater  the  difference  in  critical  visibilities  among  time  lag 
models,  and  the  more  crucial  is  the  controller  time  lag 
accuracy  for  stability. 

In  general,  if  the  guidance  law  is  slow  enough  to  allow 
sufficient  time  for  the  controller  to  react,  stability  of 
straight  line  motion  is  guaranteed  [Ref .  1] . 
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Figure  10.  Difference  of  Critical  Visibility  Distances  at 
Time  Lags  of  T  =  0.5  and  0 . 0  sec  versus  tc  for 
First,  Second  and  Third  Order  Cases. 
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Figure  11.  Difference  of  Critical  Visibility  Distances  at 
Time  Lags  of  T  =  1.0  and  0.0  sec  versus  tc  for 
First,  Second  and  Third  Order  Cases. 
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Figure  12.  Difference  of  Critical  Visibility  Distances  at 
Time  Lags  of  T  =  1.5  and  0.0  sec  versus  tc  for 
First,  Second  and  Third  Order  Cases. 
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Figure  13.  Difference  of  Critical  Visibility  Distances  at 
Time  Lags  of  T  =  2.0  and  0.0  sec  versus  tc  for 
First,  Second  and  Third  Order  Cases. 
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G.   NORMALIZATION  OF  RESULTS 

In  an  attempt  to  characterize  the  stability  analysis 
results  for  all  values  of  time  lag,  we  proceed  as  follows:  If 
the  vehicle  side  slip  velocity  v  is  very  small  so  that  it  can 
be  neglected,  the  equations  of  motion  (2.7a),  (2.7b),  (2.7c), 
and  (2.7d),  take  the  form 


i|r=r 


r=a22ur+b2u26 


y=ui|r 


since   v   =    0 .       The    control    law   is 

fi=ic1(i|r-i(;c)  +k3r, 

where  the  first  order  approximation  for  the  commanded  heading 
angle  is 

Vc    d 

Based  on  the  above  equations,  we  can  get  the  characteristic 
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equation  of  the  system  as 


s3-(a22u+b2k3u2)s2-(b2k1u2-b2k1u3T±)s-b2k1u3^  =  0 


Applying    Routh's    criterion    to    this    cubic    equation,     we    find 
that    loss    of    stability   occurs   when 


ia22u+b2k3u2)  [b2k^u2-b2k^u3T^)    =  -b2klU^, 


from  which  we   can   get    the    critical   visibility   distance 


d=- ]— —    +  Tu  (3.27) 


Equation  (3.27)  can  be  written  as 


d?   d°  =  Ui  (3.28) 


where  d0  is  the  critical  visibility  distance  for  T  =  0,  and  dT 
the  critical  visibility  distance  for  a  nonzero  T.  In 
dimensionless  quantities,  equation  (3.28)  is 
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d?  d°    =1  (3  .29) 


where  u  =  2.0  ft/sec  is  the  vehicle  speed,  and  C  =  7.3  ft  is 
the  vehicle  length. 

A  comparison  among  the  approximate  expression  (3.29)  and 
the  first,  second,  and  third  order  approximations  is  shown  in 
figures  14,  15  and  16.  The  agreement  is  better  for  the  first 
order  approximation,  as  expected,  and  also  for  the  higher 
controller  time  constant  tc .  This  is  because  a  high  tc  results 
in  soft  vehicle  response  with  negligible  amounts  of  side  slip. 


43 


2.5- 


1.5  - 


1  - 


0.5  - 


1    \\\\  '            !            ! 

■™i 

- 

, : J \>A              L 

1 

i 

\  ^ 

I 

i 

! 

\ 

1               !               \vx 

i i    \\i     i 

! 

!           0.! 

1 

i 

hi 

5  l! 

2.0 

! 

\\\ 

j 
i 

' 

1 " 

! 

1                     i 

i 

i 

i 

1 

! 

! 

i 

0.92   0.94   0.96   0.98 


1.02 


1.04    1.06    1.08 


(dT-dO)/T 


1.1 


Figure  14.  Normalized  Critical  Visibility  versus  tc  for 
First  Order  Approximation  For  Time  Lags 
T  =  0.5,  1.0,  1.5,  2.0  sees. 
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Figure  15.  Normalized  Critical  Visibility  versus  tc  for 
Second  Order  Approximation  For  Time  Lags 
T  =  0.5,  1.0,  1.5,  2.0  sees. 
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Figure  16.  Normalized  Critical  Visibility  versus  tc  for 
Third  Order  Approximation  For  Time  Lags 
T  =  0.5,  1.0,  1.5,  2.0  sees. 
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IV.   FREQUENCY  RESPONSE  ANALYSIS 

A.  INTRODUCTION 

The  previous  analysis  using  Taylor  expansions  of  y(t-T) 
breaks  down  for  approximations  beyond  third  order.  The  reason 
for  this  is  that  for  higher  than  third  order  expansions,  the 
matrix  B  in  the  generalized  eigenvalue  problem  (3.23)  becomes 
singular.  Therefore,  it  was  felt  that  a  different  technique 
should  be  sought  in  order  to  obtain  an  exact  computation  of 
the  stability  curves  and  also  to  check  the  validity  of  our 
calculations.  In  this  chapter,  we  present  a  technique  based 
on  frequency  response  methods,  which  utilizes  the  Nyquist 
criterion  for  stability.  The  resulting  equation  is  solved 
using  Newton's  iteration. 

B.  NYQUIST  CRITERION 

From  our  previously  developed  model,   the  linearized 
horizontal  equations  of  motion  are 

i|r=r  (2.7a) 

v~=a11uv+a12ur+jb1u25  (2.7b) 
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r=a12uv+a22ur  +b2u2b  (2.7c) 


y=wty  +  v  (2  .  7d) 


at  any  nominal  u.   The  control  law  is 


6=kx(y-tyc)+k2v+k3r  (2.10) 


where  the  commanded  heading  angle  is 


tc 


The  Laplace  transform  of  equation  (3.9)  is 


,-yJ^zIL  o.9) 


i|/c  =  -±e-Ts  (4.1) 


and  the  resultant  rudder  control  law  (2.10)  becomes 


b=k^+k2v+kir+-±ye-T9  (4.2) 


Following  some  algebra,  the  resultant  linearized  equations  of 
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motion  (2.7a),  (2.7b),  (2.7c),  (2.7d)  become 


t|;=r 


(4.3a) 


v=  (251u2Jc1-i?1u3^D  i|;+  (a^u+b^k^b^^T)  v+  (a^u+b^k,)  r+b^^ye-™     (4.3b) 


r=(b2u2k1-b2u1^T)  v|»+  (a21u+b2u2k2-b2u2^T)  v+  (a22u+b2u2k3)  r+b2u2^ye-Ts      (  4  .  3c 


y=in|r  +  V 


4.3d) 


Equations  (4.3a),  (4.3b),  (4.3c),  (4.3d)  reduce  to  the  state 
space  form 


x-Ax 


3.16 


The  matrix  form  becomes 


(Jb^u'leJ     (a11u+Jb1u2/c2)     (a12u+i)1u2ic3)     (i^u2  — £e 


2  ^i  ~-rs< 


(Jb^lq)     (a21u+£2u2£2)     (a22u+Jb2u2Jc3)     (£2u2-ie-rs) 
u  1  0 


d 

0 


* 


4.4 
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The   characteristic    equation   of    the    form      [A-Is]    =    0      becomes 
after   a    considerable   amount    of    algebra 

s*+A2s3+A2s2+A1s+{B2s2+B1s+B0)De-Ts  =  0  (4.5) 

where 

A2  =  -(axl+a22)u  -  [brk2+bzkz)  u2 

Ai  =  "  (312^2-^22^1)  ui^De-^-ia^-a^b^u3^ 
B2  =  -bxu2kx 


B1  =  b2u3kt 


Bo  =   (aiA^A)"4*! 


and 


--^ 
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The  characteristic  equation  of  the  system  is  written  as 


1  +  K  G(s)    =0  (4.6 


where 


G(s)  .  (B2s^Bls,B0)e-^  (4_7 

s{s2+A3s2+A2s+A1) 


is  the  transfer  function,  and  we  have  denoted 

K=D  (4.8) 

With  s  =  j(0,  the  phase  angle  <J)  is  represented  by 


<|>  =  ZG(j'to)  (4.9 


and  it  follows  that 


4>  =  /.(e-Tj(*)    -  Z(ja>)    +  Z(-B2to2  +B1jg>+.B0)    -  K-j^-A^^+A^fD+A^ 


<t>  =  -or  -  —  +  tan_1( x- )    -  tan_1( *— )    (4.10) 

2  B0-B2(D2  -A3(ji2+A1 
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For  stability,  the  Nyquist  criterion  states  (Figure  17)  that 
at  the  solution  of  the  equation 

(J)  (CD)  =  -71  (4.11) 

which  is  phase  cross  over  frequency  0)  =  0^,  the  gain  margin 
must  be  equal  to  1 

|iCG(jo>1)  |=1  (4.12) 

where   K  =  D  =  1/d.   It  follows  that  equation  (4.11)  becomes 


-c^r-iUtairM  B^    )  -tan'M  "";^Bl )  =  -*   (4.13 
2        B0-B2O)i         -A2(x>l+A1 


Rearranging,  we  get 


tan_1( ^-^)-tan_1( — ?— i)    =  -7t  +  — -t-G>r 

S0-S2a>i  -A3a>J+Z1  2 


Taking  the  tangent  of  both  sides  and  using  the  identity 


.   ,    v   tan  (x) -tan  (y) 

tan  (x-y)  = ;  ; -^-L- 

l+tan(x) tan(y) 
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GMju) 


Figure  17.  Example  of  the  Nyquist  Stability  Criterion  For 
Three  Different  Values  of  Gain  [Ref.  7]. 
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and  rearranging  following  some  algebra,  the  result  becomes 


2\  /   „3 


B1(ji1(AL-A3ta1)  -(S0-B2«i)  (-co1+A2o)1)         i 


(B0-B2g>2x)  {-A^l+Ai) +B1<*1(-«>l+A2t*1)  tan(WlD  (4.14 


From  equation  (4.12),  it  follows  that 


1  -B2v>l+Bju>i+B0\  = 


"    &)1|-ja)i-A3Ci)i+A2ja)1+A1 


Solving  for  d,  we  get 


^  _   v  v-"o    ■i^2wi/      "i  wi (4.15) 

al}J  (Ax-A3u>l)  2+  (A2u1-<j)1)  2 


With   a    time    lag   T   =    0,    we   get    from  equation    (4.14 


(fl0-B2&>i)  (^0)1+^)    +  fl1d>1(-(i>J+A2w1)    =  0 


Rearranging,  we  get 


{B2A3-Bx)a>\   +    (B^-B^-fl^)  a)2   +  B0^  =  0  (4.16 
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Equation  (4.16)  can  be  solved  exactly  for  G^,  and  then  d 
can  be  computed  from  (4.15) .  This  zero  time  lag  calculation 
was  performed  in  order  to  validate  the  results  of  d  versus  tc 
for  the  zero  time  lag  solution  of  equation  (3.6) .  The  results 
using  the  two  different  techniques  were  found  to  be  identical. 

C.   ALGORITHM  DESCRIPTION 

Introducing  the  terms 

a2  =  £iArfl2Ai  -5cA 


and 


a3  =  BQAX 


Pi   =   ~B2 


P3  =  ¥rflo^ 


p2  =  Bq+B^-B^ 
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Then   equation    (4.12)    becomes 


P1a)5  +  p2g)^P3a)    _  _  i  (41?) 

a1a)4+a2a)2+a3  tan(wr) 


Rearranging,  we  get 


(P1ci>s  +  P2g>3  +  P3&>)  tanfcor)    +  a1o)4+a2o)2  +  a3  =  0         (4.18) 


Equation    (4.18)    is    now   in   the    form 


f(W)    =  0  (4.19) 


Following  Newton's  iteration  for  CO,  we  have 


From  equation  (4.20), 


f(G>)    =    (p1co5+p2G)3  +  P3a))  sin(wT) 
+   (a1o)4+a2to2+a3)  cos  (cor) 


(4.21) 
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and 


f'(&>)    =   (5p1a)4  +  3P2G)2  +  P3)sin(a)r) 

+   (61(i)5  +  67a)2  +  B,a))  TcosiuT) 

(4.22) 

+  (4a1(i)3+2a2o))  cos  (ur) 
-  (a10><l+a2G)2+O3)  rsin(a>r) 

Fortran  programming  and  MATLAB  were  used  in  order  to 
develop  the  stability  characteristics  for  the  Newton's 
iteration  method.  Appendix  D  presents  the  program  used  to 
develop  the  stability  characteristics  for  this  case. 

Figure  18  shows  the  resulting  stability  curves  for  the 
Newton's  iteration  method.  The  results  of  the  critical 
visibility  d  versus  the  controller  time  constant  tc  for  time 
lags  of  0.0,  0.5,  1.0,  1.5,  and  2.0  sees  are  shown. 

As  with  the  first,  second  and  third  order  approximation 
cases,  the  higher  values  of  tc  require  higher  lookahead 
distances  d  for  path  accuracy,  high  values  of  d  correspond  to 
a  slow  guidance  law  with  loss  in  speed  of  response  and  path 
accuracy,  and  the  stability  curves  shift  to  the  right  with 
increasing  values  of  time  lag.  However,  since  the  Newton's 
iteration  method  presents  the  exact  solution  of  the  frequency 
response  to  equation  (4.19),  the  resulting  stability  curves 
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Figure  18.  Newton's  Iteration  Method:  Critical  Visibility 
Distance  d  versus  tc  for  Time  Lags  T  =  0.0, 
0.5,  1.0,  2.0  sees . 
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represent  the  precise  stability  characteristics  for  these  time 
lags  . 

D.   RESULTS  AND  DISCUSSION 

Results  of  the  critical  visibility  d  versus  the  controller 
time  constant  tc  for  the  Newton's  iteration  method  has  been 
demonstrated.  Figure  18  presented  the  precise  stability 
curves  for  time  lags  of  0.0,  0.5,  1.0,  1.5  and  2.0  sees. 

This  figure,  as  well  as  those  for  the  first,  second  and 
third  order  approximation  cases,  demonstrates  that  the 
stability  of  the  control  scheme  decreases  with  increase  in 
time  lag  and  for  the  same  time  constant,  critical  visibility 
increases  with  increase  in  time  lag. 

Figures  19,  20,  21,  22  and  23  present  a  comparison  of 
generated  stability  curves  among  the  Newton's  iteration  method 
and  those  for  the  first,  second  and  third  order  approximation 
cases  for  time  lags  of  0.0,  0.5,  1.0,  1.5  and  2.0  sees.  It 
can  be  seen  that  there  is  barely  a  slight  difference  in  the 
stability  curves  for  the  Newton's  iteration  method  case  and 
the  third  order  approximation  case.  This  indicates  that  a 
third  order  approximation  for  time  lag  presents  the  best  model 
for  accounting  for  a  positional  information  time  lag  in  the 
vehicle  control  law. 
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Figure  19.  Critical  Visibility  Distance  d  versus  tc  for 
Newton's  Iteration  Method  and  First,  Second  and 
Third  Order  Approximations;  T  =  0.0  sec. 
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Figure  20.  Critical  Visibility  Distance  d  versus  tc  for 
Newton's  Iteration  Method  and  First,  Second  and 
Third  Order  Approximations;  T  =  0.5  sec. 


61 


Figure  21.  Critical  Visibility  Distance  d  versus  tc  for 
Newton's  Iteration  Method  and  First,  Second  and 
Third  Order  Approximations;  T  =  1.0  sec. 
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Figure  22.  Critical  Visibility  Distance  d  versus  tc  for 
Newton's  Iteration  Method  and  First,  Second  and 
Third  Order  Approximations;  T  =  1.5  sec. 
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d  vs  tc  for  1st,  2nd,  3rd  Ord  &  nefrtons  method  time  lag  T=2.0 


Figure  23.  Critical  Visibility  Distance  d  versus  tc  for 
Newton's  Iteration  Method  and  First,  Second  and 
Third  Order  Approximations;  T  =  2.0  sec. 
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CONCLUSIONS  AND  RECOMMENDATIONS 

A.   CONCLUSIONS 

A  methodology  for  considering  positional  information  time 
lags  in  the  control  law  for  vehicle  maneuvering  in  the 
horizontal  plane  has  been  presented.  The  relationship  of  the 
critical  visibility  versus  the  controller  time  constant  for 
both  zero  and  non-zero  time  lags  has  been  established. 

Time  lags  have  been  approximated  by  first,  second  and 
third  order  Taylor  expansions  in  the  commanded  vehicle  heading 
in  the  control  law.  A  comparison  of  the  resulting  stability 
curves  has  demonstrated  that  the  third  order  approximation 
presents  the  greatest  stability  and  the  least  critical 
visibility  for  a  given  time  lag.  It  was  found  that  the  first 
order  approximation  was  the  most  conservative  for  controller 
design  purposes  since  it  predicted  the  highest  required 
critical  visibility  distance.  In  addition,  the  differences 
among  the  three  approximations  became  more  pronounced  with 
increasing  time  lag  and  decreasing  controller  time  constant 
(tighter  control). 

Further  analysis  was  conducted  to  characterize  the 
stability  results  for  any  value  of  time  lag  using  a 
normalization  procedure  and  the  first  order  approximation  for 
time  lag  in  the  commanded  vehicle  heading  in  the  control  law. 
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A  comparison  of  the  normalized  stability  curves  for  each  time 
lag  showed  that  the  greatest  differences  occur  with  increasing 
order  and  decreasing  controller  time  constant.  This  is  due  to 
a  tight  (quick)  vehicle  response  with  significant  side  slip  at 
lower  time  constants. 

Frequency  response  techniques  based  on  the  Nyquist 
stability  criterion  were  used  to  produce  the  exact  stability 
curves  for  comparison  with  the  previously  generated  Taylor 
expansion  stability  curves.  These  results  virtually  matched 
those  obtained  for  the  third  order  approximation  time  lag  case 
and  validated  our  calculations  and  stability  curves  generated 
for  the  case  of  zero  time  lag.  These  results  indicated  that 
the  third  order  approximation  best  represents  a  positional 
information  time  lag  in  the  control  law,  however  for  design 
purposes,  the  first  order  approximation  provides  the  most 
conservative  approach. 

B .   RECOMMENDATIONS 

Some  recommendations  for  further  research  are  as  follows: 

*  Experimental  verification  using  the  NPS  AUV  II. 

*  Consideration  of  time  lag  in  the  simulation  of  an  Inertial 
Navigation  System  required  for  positional  updates. 

*  Analysis  of  positional  information  time  lags  in  the  motion 
stability  in  the  vertical  plane,  along  curved  paths,  or 
with  respect  to  other  guidance  schemes. 


66 


APPENDIX    A 


PROGRAM    THESISl.FOR    (Time    Delay-lst    Order    Approx) 
BIFURCATION    ANALYSIS 

IMPLICIT    DOUBLE    PRECISION     (A-H,0-Z) 

DOUBLE    PRECISION    Kl , K2 , K3 , L , NR , NV , NDRS , NDRB , I Z , MASS , 
&  NRDOT,NVDOT 

DIMENSION    A(4,4),FVl(4),IVl(4),ZZZ(4,4),WR<4),WI(4) 

OPEN  (11, FILE='BIFl .RES' ,STATUS='NEW ) 
OPEN  (12, FILE='BIF2.RES' ,STATUS='NEW ) 
OPEN    (13,FILE='BIF3 .RES' , STATUS= ' NEW ) 

Vehicle  Parameters: 
WEIGHT=435.0 


IZ 
L 

RHO 
G 

XG 
MASS 


=  45.0 
=  7.3 
=  1.94 
=  32.2 
=0.0104 
=WEIGHT/G 


YRDOT  =-0 

YVDOT  =-0 

YR 

YV 

YDRS 

YDRB 

NRDOT 

NVDOT 

NR 

NV 

NDRS 

NDRB 


=  +  0 

—  0 
=  +  0 
=  +  0 
=  -0 
=  -0 
=  -0 

— o 

=-0 

-  +  0 


00178*0. 5*RHO*L**4 
03430*0. 5*RHO*L**3 
01187*0. 5*RHO*L**3 
03896*0. 5*RHO*L**2 
02345*0. 5*RHO*L**2 
02345*0. 5*RHO*L**2 
00047*0. 5*RHO*L**5 
00178*0.  5*RHO*L**4 
01022*0. 5*RHO*L**4 
00769*0. 5*RHO*L**3 
337*0.02345*0. 5*RHO*L**3 
283*0.02345*0. 5*RHO*L**3 


WRITE  (* 

READ  ( * 

WRITE  (* 

READ  ( * 


1001 
*) 

1002 
*) 


TCMIN.TCMAX, ITC 
XDMIN.XDMAX, IXD 


XDMIN=XDMIN*L 
XDMAX=XDMAX*L 
WRITE  (*,1003) 
READ   ( * , * ) 


WRITE 
READ 


DH 

AA1 
AAl 

AA2 
AA2 

BB1 
BBl 

BB2 
BB2 


*,1100) 
*,*)  TL 

I Z-NRDOT ) * ( MASS- YVDOT ) - 

MASS*XG-YRDOT) * ( MASS *XG-NVDOT ) 

( IZ-NRDOT) *YV-(MASS*XG-YRDOT) *NV)/DH 

(  IZ-NRDOT)  M-MASS  +  YR)- 

MASS*XG-YRDOT) * ( -MASS*XG+NR) )/DH 

(MASS- YVDOT )*NV-( MASS *XG-NVDOT)*YV)/DH 

( MASS- YVDOT)* (-MASS* XG+NR)- 

MASS*XG-NVDOT) * ( -MASS+YR) )/DH 

( I Z-NRDOT )* YDRS- ( MASS *XG- YRDOT) *NDRS)/DH 

( I Z-NRDOT)* YDRB- ( MASS *XG- YRDOT)* NDRB )/DH 

(MASS- YVDOT) * NDRS- ( MASS *XG-NVDOT ) *YDRS)/DH 

(MASS- YVDOT) *NDRB- ( MASS *XG-NVDOT )* YDRB ) /DH 


WRITE  (*,1004) 
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READ   (*,*)      RATIO 

BB1=BB11+RATI0*BB12 
BB2=BB21+RATIO*BB2  2 

EPS   =l.D-5 
ILMAX=1500 

DO  1  1=1 , ITC 

WRITE  (*,2001)  I,ITC 

TC=TCMIN+( 1-1 )*(TCMAX-TCMIN)/( ITC-1 ) 

TCD=TC*L/U 

POLEC=l . 0/TCD 

ADl-3 . 0*POLEC 

AD2=3 . 0*POLEC**2 

AD3=POLEC**3 

A1=BB1*U*U 

B1=BB2*U*U 

Cl=-ADl-(AAll+AA2  2 ) *U 

A2=(BB1*AA22-BB2*AA12)*U**3 

B2=(BB2*AAll-BBl*AA21 )*U**3 

Kl=AD3/( (BB2*AA11-BB1*AA21)*U**3) 

C2=AD2-( AA11*AA22-AA12*AA21 )*U**2+BB2*U*U*Kl 

K2=(C1*B2-C2*B1 )/( Al *B2-A2 *Bl ) 

K3=(C2*A1-C1*A2)/(A1*B2-A2*B1) 

DO  2  J=l , IXD 

XD=XDMIN+( J-l ) * (XDMAX-XDMIN)/( IXD-1 ) 

CALL  LINEAR (Kl , K2 , K3 , U , XD, AAl 1 , AAl 2 , AA21 , AA2  2 , BBl ,BB2,A,TL) 
CALL  RG(4,4,A,WR,WI,0,ZZZ,IV1,FV1, I  ERR) 
CALL  DSTABL(DEOS,WR,WI , FREQ) 

IF  ( J .GT. 1 )  GO  TO  10 
DEOSOO=DEOS 
XDOO   =XD 
LL  =  0 
GO  TO  2 
10      DEOSNN-DEOS 
XDNN   =XD 
PR=DEOSNN*DEOSOO 
IF  ( PR.GT.0.D0)  GO  TO  3 
LL=LL+1 

IF  (LL.GT.3)  STOP  1000 
IL  =  0 

XDO=XDOO 
XDN=XDNN 
DEOSO=DEOSOO 
DEOSN-DEOSNN 
6       XDL=XDO 
XDR=XDN 
DEOSL=DEOSO 
DEOSR=DEOSN 
XD=(XDL+XDR)/2 .DO 

CALL  LI  NEAR ( Kl,K2,K3,U,XD,AAll ,AAl2, AA21 ,AA2  2,BBl ,BB2,A,TL) 
CALL  RG(4(4,A/WR,WI,0,ZZZ,IV1,FV1, I  ERR) 
CALL  DSTABL( DEOS,WR,WI , FREQ) 

DEOSM=DEOS 
XDM=XD 
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PRL=DEOSL*DEOSM 
PRR=DEOSR*DEOSM 
IF  ( PRL.GT. 0 .DO)  GO  TO  5 
XDO=XDL 
XDN=XDM 
DEOSO=DEOSL 
DEOSN=DEOSM 
IL=IL+1 

IF  ( IL.GT. ILMAX)  STOP  3100 
DIF=DABS(XDL-XDM) 
IF  (DIF.GT.EPS)  GO  TO  6 
XD=XDM 
GO  TO  4 

IF  (PRR.GT.0.D0)  STOP  3200 
XDO=XDM 
XDN=XDR 
DEOSO=DEOSM 
DEOSN=DEOSR 
IL=IL+1 

IF  ( IL.GT. ILMAX)  STOP  3100 
DIF=DABS(XDM-XDR) 
IF  (DIF.GT.EPS)  GO  TO  6 
XD=XDM 
LLL=10+LL 

WRITE  (LLL,*)  XD/L,TC 
XDOO=XDNN 
DEOSOO=DEOSNN 
CONTINUE 
CONTINUE 


1001  FORMAT  ( 

1002  FORMAT  ( 
100  3  FORMAT  ( 
1004  FORMAT  ( 
1100  FORMAT  ( 
2001  FORMAT  (215) 

END 


ENTER  MIN,  MAX,  AND  INCREMENTS  OF  TC ' ) 

ENTER  MIN,  MAX,  AND  INCREMENTS  OF  XD'  ) 

ENTER  U'  ) 

ENTER  BOW/STERN  RUDDER  RATIO' ) 

ENTER  TIME  LAG  TL' ) 


SUBROUTINE  DSTABL ( DEOS , WR , WI , OMEGA) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  WR( 4 ) ,WI( 4) 
DEOS--1 . 0D+20 
DO  1  1=1 , 4 

IF  (WR( I ) .LT.DEOS)  GO  TO  1 

DEOS=WR( I ) 

IJ=I 
1  CONTINUE 

OMEGA-WI ( IJ) 
OMEGA=DABS ( OMEGA ) 
RETURN 
END 

SUBROUTINE  LINEAR(Kl , K2 , K3 , U , XD , AAl 1 , AAl 2 , AA21 , AA2  2 , BBl , BB2 , A, TL ) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  Kl,K2,K3 

DIMENSION  A( 4,4) 

A(  1,  D-0.0D0 

A(1,2)=0.0D0 

A(1,3)=1.0D0 

A(1,4)-0.0D0 
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A(2,l)»  BB1*U*U*K1*(XD-U*TL)/XD 

A(  2,2)=U*(BB1*U*(K2*XD-K1*TL)+AA11*XD)/XD 

A(  2,  3  )=AAl2*U  +  BBl*U*U*K3 

A(2,4)=  BB1*U*U*K1/XD 

A(3,l)=»  BB2*U*U*Kl*(XD-U*TL)/XD 

A( 3, 2)=U*(BB2*U*(K2*XD-K1*TL)+AA21*XD)/XD 

A(  3,  3 )=AA2  2*U+BB2*U*U*K3 

A(3,4)»  BB2*U*U*Kl/XD 

A(4,1)=U 

A(  4,2)=1 .ODO 

A(4,  3)=O.ODO 

A( 4 , 4 )=0 . ODO 

RETURN 

END 
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APPENDIX    B 


PROGRAM   THESIS2.FOR    (Time    Delay-2nd    Order    Approx) 
BIFURCATION    ANALYSIS 

IMPLICIT    DOUBLE    PRECISION    (A-H,0-Z) 

DOUBLE    PRECISION    Kl , K2 , K3 , L , NR , NV , NDRS , NDRB , I Z , MASS , 
&  NRDOT,NVDOT 

DIMENSION    A(4,4) , FVl ( 4 ) ,IVl(4),ZZZ(4,4),WR(4),WI(4) 


OPEN  (ll,FILE='BIFl .RES' , STATUS= ' NEW ) 
OPEN  ( 12,FILE='BIF2 .RES' ,  STATUS='NEW ) 
OPEN    (13,FILE='BIF3.RES' ,  STATUS- 'NEW ) 

Vehicle  Parameters: 
WEIGHT-435. 0 


IZ 

L 

RHO 

G 

XG 

MASS 


=  45.0 
=  7.3 
=  1.94 
=  32.2 
=0.0104 
=WEIGHT/G 


YRDOT  =-0. 00178*0. 5*RHO*L**4 
YVDOT  =-0. 03430*0. 5*RHO*L**3 
YR  =+0.01187*0. 5*RHO*L**3 
YV  =-0.03896*0. 5*RHO*L**2 
YDRS  =+0.02345*0. 5*RHO*L**2 
YDRB  =+0.02345*0. 5*RHO*L**2 
NRDOT  =-0.00047*0. 5*RHO*L**5 
NVDOT  --0. 00178*0. 5*RHO*L**4 
NR  --0.01022*0. 5*RHO*L**4 
NV  =-0. 00769*0. 5*RHO*L**3 
NDRS  =-0.  3  37*0.02  345*0.5*RHO*L**3 
NDRB   =+0.28  3*0.02  345*0.5*RHO*L**3 


WRITE  (*,1001) 
READ   (*,*) 
WRITE  (*,1002) 
READ   (*,*) 
XDMIN=XDMIN*L 
XDMAX=XDMAX*L 
WRITE  (*,1003) 
READ   (*,*) 


TCMIN.TCMAX, ITC 
XDMIN,XDMAX, IXD 


WRITE 

[*,1100) 

READ 

*,*)  TL 

DH   = 

IZ-NRDOT)*(MASS 

MASS*XG-YRDOT) * 

AA11- 

( IZ-NRDOT)*YV-( 

AA12- 

(  IZ-NRDOT) *(-MA 

MASS*XG-YRDOT)* 

AA21- 

(MASS-YVDOT) *NV 

AA22- 

(MASS-YVDOT)*(- 

MASS*XG-NVDOT)* 

BB11- 

( IZ-NRDOT) *YDRS 

BB12- 

( IZ-NRDOT) *YDRB 

BB21  = 

(MASS-YVDOT) *ND 

BB22  = 

(MASS-YVDOT) *ND 

-YVDOT)- 

(MASS*XG-NVDOT) 

MASS*XG-YRDOT) *NV)/DH 

SS+YR)- 

(-MASS*XG+NR) )/DH 

-(MASS*XG-NVDOT)*YV)/DH 

MASS*XG+NR)- 

(-MASS+YR) )/DH 

- ( MASS *XG- YRDOT ) *NDRS )/DH 

-( MASS *XG- YRDOT)* NDRB )/DH 

RS-(MASS*XG-NVDOT) *YDRS)/DH 

RB-(MASS*XG-NVDOT) *YDRB)/DH 


WRITE  (*,1004) 
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READ   (*,*)      RATIO 

BB1=BB11+RATI0*BB12 
BB2=BB21+RATIO*BB2  2 

EPS   =l.D-5 
ILMAX=1500 

DO  1  1=1 , ITC 

WRITE  (*,2001)  I,ITC 

TC=TCMIN+( 1-1 ) *(TCMAX-TCMIN)/( ITC-1 ) 

TCD=TC*L/U 

POLEC-1 . 0/TCD 

ADl-3 . 0*POLEC 

AD2=3.0*POLEC**2 

AD3=POLEC**3 

Al=BBl*U*U 

Bl=BB2*U*U 

Cl»-ADl-( AAll+AA2  2)*U 

A2=(BB1*AA22-BB2*AA12)*U**3 

B2=(BB2*AAll-BBl*AA21 )*U**3 

Kl=AD3/( (BB2*AA11-BB1*AA21 ) *U**3) 

C2-AD2-( AA11*AA2  2-AA12*AA21 ) *U* * 2+BB2*U*U*Kl 

K2=(C1*B2-C2*B1 )/( Al *B2-A2 *Bl  ) 

K3=(C2*A1-C1*A2)/(A1*B2-A2*B1) 

DO  2  J=l , IXD 

XD»XDMIN+( J-l )*(XDMAX-XDMIN)/( IXD-1 ) 

CALL  LINEAR( Kl ,K2,K3,U,XD,AAll ,AA12, AA21 ,AA2  2,BBl ,BB2,A,TL) 
CALL  RG(4,4,A,WR,WI,0,ZZZ,IVl,FVl, I  ERR) 
CALL  DSTABL( DEOS,WR,WI , FREQ) 

IF  ( J.GT. 1 )  GO  TO  10 
DEOSOO=DEOS 
XDOO   =XD 
LL  =  0 
GO  TO  2 
10       DEOSNN=DEOS 
XDNN   =XD 
PR=DEOSNN*DEOSOO 
IF  ( PR.GT.0.D0)  GO  TO  3 
LL=LL+1 

IF  (LL.GT.3)  STOP  1000 
IL  =  0 

XDO=XDOO 
XDN=XDNN 
DEOSO=DEOSOO 
DEOSN=DEOSNN 
6      XDL=XDO 
XDR=XDN 
DEOSL-DEOSO 
DEOSR=DEOSN 
XD=(XDL+XDR)/2 .DO 

CALL  LI  NEAR ( Kl , K2 , K3 , U , XD , AAl 1 , AAl 2 , AA2 1 , AA2 2 , BBl , BB2 , A , TL ) 
CALL  RG( 4, 4,A,WR,WI ,0,ZZZ, IVl.FVl, I  ERR) 
CALL  DSTABL(DEOS,WR,WI , FREQ) 

DEOSM=DEOS 
XDM=XD 
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PRL-DEOSL*DEOSM 

PRR=DEOSR*DEOSM 

IF  ( PRL.GT.O.DO)  GO  TO  5 

XDO=XDL 

XDN=XDM 

DEOSO=DEOSL 

DEOSN=DEOSM 

IL=IL+1 

IF  ( IL.GT. ILMAX)  STOP  3100 

DIF=DABS(XDL-XDM) 

IF  (DIF.GT.EPS)  GO  TO  6 


XD=XDM 

GO  TO  4 

5 

IF  (PRR.GT.O.DO)  STOP 

XDO=XDM 

XDN=XDR 

DEOSO=DEOSM 

DEOSN-DEOSR 

IL=IL+1 

3200 

IF  ( IL.GT. ILMAX)  STOP 

3100 

DIF=DABS(XDM-XDR) 

IF  (DIF.GT.EPS)  GO  TO 

6 

XD=XDM 

4 

LLL=10+LL 

WRITE  (LLL,*)  XD/L,TC 

3 

XDOO=XDNN 
DEOSOO=DEOSNN 

2 

CONTINUE 

1 

CONTINUE 

1001 

FORMAT  ( '  ENTER  MIN,  MAX, 

AND  INCREMENTS 

OF  TC 

1002 

FORMAT  ('  ENTER  MIN,  MAX, 

AND  INCREMENTS 

OF  XD 

1003 

FORMAT  ( '  ENTER  U' ) 

1004 

FORMAT  ( '  ENTER  BOW/STERN 

RUDDER  RATIO' ) 

1100 

FORMAT  ( '  ENTER  TIME  LAG  TL' ) 

2001 

FORMAT  (215) 
END 

SUBROUTINE  DSTABL ( DEOS , WR , WI , OMEGA) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  WR( 4) ,WI( 4 ) 
DEOS=-l .OD+20 
DO  1  1=1,4 

IF  (WR( I ) .LT.DEOS)  GO  TO  1 

DEOS=WR( I ) 

IJ=I 
CONTINUE 
OMEGA=WI( IJ) 
OMEGA=DABS ( OMEGA ) 
RETURN 
END 


SUBROUTINE  LI NEAR ( Kl , K2 , K3 , U , XD, AAl 1 , AAl 2 , AA21 , AA2  2 , BBl 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  Kl,K2,K3 

DIMENSION  A( 4, 4) 

A(  1,1)=0.0D0 

A(1,2)=0.0D0 

A(1,3)=1.0D0 

A(1,4)=0.0D0 


BB2 ,A,TL) 
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A  (2,  1  )-(BBl*U*U*Kl-BBl*U*U*U*Kl*TL/XD)/ 
&  ( l-BBl*U*U*Kl*TL*TL/( 2.D0*XD) ) 

A( 2, 2)=(AAll*U+BBl*U*U*K2-BBl*U*U*Kl*TL/XD)/ 
&  ( l-BBl*U*U*Kl*TL*TL/( 2 .DO*XD) ) 

A( 2, 3)=(AAl2*U+BBl*U*U*K3+BBl*U*U*U*Kl*TL*TL/( 2.D0*XD) )/ 
&  ( l-BBl*U*U*Kl*TL*TL/( 2.D0*XD) ) 

A(  2,  4  )  =  (BB1*U*U*K1/XD)/ 
&  ( l-BBl*U*U*Kl*TL*TL/( 2 .DO*XD) ) 

A(  3 , 1  )=(BB2*U*U*Kl-BB2*U*U*U*Kl*TL/XD)+ 
&  (BB2*U*U*Kl*TL*TL/( 2 .DO*XD) ) * 

&  (BBl*U*U*Kl-BBl*U*U*U*Kl*TL/XD)/ 

&  ( l-BBl*U*U*Kl*TL*TL/( 2 .DO*XD) ) 

A( 3, 2)=(AA21*U+BB2*U*U*K2-BB2*U*U*K1*TL/XD)+ 
&  (BB2*U*U*Kl*TL*TL/( 2.D0*XD) )* 

&  (AAll*U+BBl*U*U*K2-BBl*U*U*Kl*TL/XD)/ 

&  ( l-BBl*U*U*Kl*TL*TL/( 2.D0*XD) ) 

A( 3, 3 )=( AA22*U+BB2*U*U*K3+BB2*U*U*U*Kl*TL*TL/(  2 .DO*XD)  )  + 
&  (BB2*U*U*Kl*TL*TL/( 2 .DO*XD) )* 

&  ( AAl2*U+BBl*U*U*K3+BBl*U*U*U*Kl*TL*TL/( 2 .DO*XD) )/ 

&  ( l-BBl*U*U*Kl*TL*TL/( 2 .DO*XD) ) 

A( 3, 4 )=(BB2*U*U*Kl/XD)+( BB2 *U*U*Kl *TL*TL/( 2 .DO*XD) )* 
&  (BB1*U*U*K1/XD)/ 

&  ( l-BBl*U*U*Kl*TL*TL/(2.D0*XD) ) 

A(4,l)-U 

A(4,2)«l .ODO 

A(  4, 3)=O.ODO 

A(4,4)=0.0D0 

RETURN 

END 
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APPENDIX  C 


PROGRAM  THESIS3.FOR  (Time  Delay-3rd  Order  Approx) 
BIFURCATION  ANALYSIS 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  Kl , K2 , K 3 , L , NR , NV , NDRS , NDRB , I Z , MASS , 
4,  NRDOT,NVDOT 

DIMENSION  A(5,5) ,  B(5,5) ,FVl(5) , IVl(5) , ZZZ ( 5 , 5 ) , ALFR( 5 ) , 
&  ALFI( 5) ,BETA( 5) ,WR( 5) ,WI( 5) 


OPEN  (11, FILE-'BIFl .RES' , STATUS- ' NEW' ) 
OPEN  (12, FILE='BIF2.RES' ,STATUS='NEW ) 
OPEN  (13,FILE='BIF3.RES' , STATUS- ' NEW' ) 

Vehicle  Parameters: 
WEIGHT=435.0 


IZ 

=  45.0 

L 

=  7.3 

RHO 

-1.94 

G 

-32.2 

XG 

=0.0104 

MASS 

=WEIGHT/G 

YRDOT 

=-0.00178*0.5*RHO*L**4 

YVDOT 

=-0.03430*0. 5*RHO*L**3 

YR 

-+0. 01187*0. 5*RHO*L** 3 

YV 

=-0.03896*0.5*RHO*L**2 

YDRS 

-+0.02345*0. 5*RHO*L**2 

YDRB 

-+0.02  34  5*0.5*RHO*L**2 

NRDOT 

=-0.00047*0. 5*RHO*L**5 

NVDOT 

--0.00178*0.5*RHO*L**4 

NR 

--0. 01022*0. 5*RHO*L**4 

NV 

--0.00769*0. 5*RHO*L**3 

NDRS 

=-0. 3  37*0.02345*0.5*RHO*L**3 

NDRB 

=+0.28  3*0.02  345*0.5*RHO*L**3 

WRITE 

(*,1001) 

READ 

(*,*)     TCMIN,TCMAX, ITC 

WRITE 

(*,1002) 

READ 

(*,*)     XDMIN,XDMAX, IXD 

XDMIN=XDMIN*L 
XDMAX=XDMAX*L 
WRITE  (*,1003) 
READ   (*,*) 


WRITE 
READ 

DH   - 

AA11- 
AA12- 

AA21- 
AA22- 

BB11- 
BB12- 
BB21- 
BB22- 


*,1100) 
*, *)  TL 

IZ-NRDOT) *( MASS- YVDOT ) - 

MASS*XG-YRDOT)*(MASS*XG-NVDOT) 

( IZ-NRDOT) *YV-(MASS*XG-YRDOT) *NV)/DH 

( IZ-NRDOT) *(-MASS+YR)- 

MASS*XG-YRDOT) * ( -MASS*XG+NR) )/DH 

(MASS-YVDOT) *NV- ( MASS *XG-NVDOT ) * YV )/DH 

(MASS- YVDOT) * ( -MASS*XG+NR ) - 

MASS*XG-NVDOT) * ( -MASS+YR) )/DH 

( I Z-NRDOT)* YDRS- ( MASS *XG- YRDOT )*NDRS)/DH 

( IZ-NRDOT) *YDRB-(MASS*XG-YRDOT) *NDRB)/DH 

(MASS-YVDOT) * NDRS- ( MASS  * XG-NVDOT ) *YDRS)/DH 

(MASS-YVDOT) *NDRB- ( MASS  * XG-NVDOT ) *YDRB)/DH 


WRITE    (*,1004) 
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c 


READ   (*,*)      RATIO 

BBl-BBll+RATIO*BBl2 
BB2-BB21+RATIO*BB22 

EPS   -l.D-5 
ILMAX=1500 

DO    1     I=1,ITC 

WRITE    (*,2001)     I,ITC 

TC=>TCMIN+( 1-1 ) *(TCMAX-TCMIN)/( ITC-1 ) 

TCD-TC*L/U 

POLEC=l .0/TCD 

ADl=3.0*POLEC 

AD2-3 . 0*POLEC**2 

AD3=POLEC**3 

Al=BBl*U*U 

Bl-BB2*U*U 

Cl=-ADl-( AA11+AA2  2) *U 

A2=(BB1*AA22-BB2*AA12)*U**3 

B2=(BB2*AAll-BBl*AA21 )*U**3 

Kl=AD3/( (BB2*AA11-BB1*AA21 )*U**3) 

C2-AD2-(AA11*AA2  2-AA12*AA21 )*U**2+BB2*U*U*Kl 

K2=(Cl*B2-C2*Bl )/( Al *B2-A2 *Bl ) 

K3=(C2*A1-C1*A2)/(A1*B2-A2*B1) 

DO  2  J=1,IXD 

XD-XDMIN+( J-l )*(XDMAX-XDMIN)/( IXD-1) 

CALL  LINEAR(K1 , K2 , K3 , U , XD , AAl 1 , AAl 2 , AA21 , AA22 , BBl , BB2 , A, B , TL ) 
CALL  RGG( 5, 5, A, B , ALFR , ALFI ,BETA,0,ZZZ, I  ERR) 
DO  11  IJE-1,5 

WR( IJE)=ALER( IJE)/BETA( IJE) 

WI( IJE)=ALFI( IJE)/BETA( IJE) 
11      CONTINUE 

CALL  DSTABL(DEOS,WR,WI , FREQ) 

IF  ( J.GT.l )  GO  TO  10 
DEOSOO=DEOS 
XDOO   =XD 
LL  =  0 
GO  TO  2 
10      DEOSNN=DEOS 
XDNN   =XD 
PR=DEOSNN*DEOSOO 
IF  (PR.GT.0.D0)  GO  TO  3 
LL-LL+1 

IF  (LL.GT.3)  STOP  1000 
IL-0 

XDO=XDOO 
XDN=XDNN 
DEOSO«=DEOSOO 
DEOSN=DEOSNN 
6       XDL-XDO 
XDR=XDN 
DEOSL=DEOSO 
DEOSR=DEOSN 
XD=>(XDL  +  XDR)/2  .DO 

CALL  LI  NEAR ( Kl,K2,K3,U,XD,AAll,AAl2/AA21,AA22,BBl,BB2,A/B,TL) 
CALL  RGG( 5,5, A, B, ALFR,ALFI , BETA, 0,ZZZ,  I  ERR) 
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DO  12  IJE- 1,5 

WR( IJE)-ALFR( IJE)/BETA( IJE) 

WI( IJE)=ALFI ( IJE)/BETA( IJE) 
12      CONTINUE 

CALL  DSTABL( DEOS,WR,WI , FREQ) 

DEOSM-DEOS 

XDM=XD 

PRL=DEOSL*DEOSM 

PRR=DEOSR*DEOSM 

IF  ( PRL.GT.O.DO)  GO  TO  5 

XDO-XDL 

XDN=XDM 

DEOSO-DEOSL 

DEOSN-DEOSM 

IL-IL+1 

IF  ( IL.GT. ILMAX)  STOP  3100 

DIF=DABS(XDL-XDM) 

IF  (DIF.GT.EPS)  GO  TO  6 

XD=XDM 

GO  TO  4 
5      IF  (PRR.GT.0.D0)  STOP  3200 

XDO=XDM 

XDN-XDR 

DEOSO=DEOSM 

DEOSN=DEOSR 

IL=IL+1 

IF  ( IL.GT. ILMAX)  STOP  3100 

DIF-DABS(XDM-XDR) 

IF  (DIF.GT.EPS)  GO  TO  6 

XD=XDM 
4      LLL=10+LL 

WRITE  (LLL,*)  XD/L,TC 
3      XDOO-XDNN 

DEOSOO-DEOSNN 
2    CONTINUE 
1  CONTINUE 

■ 

1001  FORMAT  ('  ENTER  MIN,  MAX,  AND  INCREMENTS  OF  TC) 

1002  FORMAT  ('  ENTER  MIN,  MAX,  AND  INCREMENTS  OF  XD'  ) 

1003  FORMAT  ('  ENTER  U') 

1004  FORMAT  ('  ENTER  BOW/STERN  RUDDER  RATIO') 
1100  FORMAT  ('  ENTER  TIME  LAG  TL') 

2001  FORMAT  (215) 
END 

SUBROUTINE  DSTABL ( DEOS , WR , WI , OMEGA) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  WR(5),WI(5) 
DEOS--1 .OD+20 
DO  1  1=1,5 

IF  (WR( I  )  .LT.DEOS)  GO  TO  1 

DEOS=WR( I ) 

IJ-I 
1  CONTINUE 

OMEGA=WI( IJ) 
OMEGA=DABS ( OMEGA ) 
RETURN 
END 
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SUBROUTINE  LI NEAR ( Kl , K2 , K3 , U , XD , AAl 1 , AA12 , AA21 , AA2 2 , BBl ,BB2 
A,B,TL) 
BLE  PRECISION  (A-H.O-Z) 
SION  Kl,K2,K3 
5,5),B(5,5) 


IMPLICI 

DOUBLE 

DIMENSI 


A( 

1  , 

1) 

A( 

1  , 

2) 

A( 

1  , 

3) 

A( 

1, 

4) 

A( 

1 

5) 

A( 

2 

1) 

A( 

2 

2) 

A| 

2 

3) 

A 

2 

4) 

A( 

2 

5) 

A( 

3 

1) 

A( 

3, 

2) 

A( 

3 

3) 

Al 

3 

4) 

A< 

3 

5) 

Al 

4 

1) 

A 

4 

2) 

A 

4 

3) 

A( 

4 

4) 

A( 

4 

5) 

A 

5 

1) 

A( 

5 

2) 

A( 

5 

3) 

A( 

5 

4) 

A( 

5 

5) 

B( 

1, 

1  ) 

B( 

1 

2) 

B 

1 

3) 

B 

1 

4) 

B 

1 

5) 

B 

2 

1) 

B 

2 

2) 

B 

2 

3) 

B 

2 

4) 

B 

2 

5) 

B 

3 

1) 

B 

3 

2) 

B 

3 

3) 

B 

3 

4) 

B 

3 

5) 

B 

4 

1) 

B 

4 

2) 

B 

4 

3) 

B 

4 

4) 

B 

4 

5) 

B 

5 

1) 

B 

5 

2) 

B 

5 

3) 

B 

5 

4) 

B 

5 

5) 

T  DOU 

PRECI 

ON  A( 

0.0D0 

0.0D0 

1  .0D0 

0.0D0 

0.0D0 

0.0D0 

0.0D0 

0.0D0 

0.0D0 

1.0D0 

(BBl* 

(AAll 

(AA12 

(BBl* 

(-1.D 

U 

1  .0D0 

0.0D0 

O.ODO 

0.0D0 

(BB2* 

(AA21 

(  AA22 

(BB2* 

(  BB2* 


U*U*Kl-BBl*U*U*U*Kl*TL/XD) 
*U+BB1*U*U*K2-BB1*U*U*K1*TL/XD) 
*U+BBl*U*U*K3+BBl*U*U*U*Kl*TL*TL/( 2 .DO*XD) ) 
U*U*Kl/XD) 
0+BBl*U*U*Kl*TL*TL/( 2 .D0*XD) ) 


U*U*Kl-BB2*U*U*U*Kl*TL/XD) 
*U+BB2*U*U*K2-BB2*U*U*Kl*TL/XD) 
*U+BB2*U*U*K3+BB2*U*U*U*Kl*TL*TL/( 2 
U*U*K1/XD) 
U*U*K1*TL*TL/(2.D0*XD) ) 


DO*XD)  ) 


=1 .ODO 

-O.ODO 

-O.ODO 

=0.0D0 

=0.0D0 

-O.ODO 

-1.0D0 

-O.ODO 

-O.ODO 

-O.ODO 

-O.ODO 

-O.ODO 

=(BBl*U*U*U*Kl*TL*TL*TL/(6 . DO*XD) ) 

-O.ODO 

-( BBl*U*U*Kl*TL*TL*TL/(6 . D0*XD) ) 

-O.ODO 

-O.ODO 

-O.ODO 

-1 .ODO 

-O.ODO 

-O.ODO 

-O.ODO 

-( 1 .0D0+BB2*U*U*U*K1*TL*TL*TL/(6.D0*XD) ) 

-O.ODO 

-(BB2*U*U*K1*TL*TL*TL/(6.D0*XD)) 


RETURN 
END 
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APPENDIX    D 


PROGRAM   THESISN.FOR    (Time    Delay    -    Newton's    Iteration    Method) 

IMPLICIT    DOUBLE    PRECISION     (A-H,0-Z) 

DOUBLE    PRECISION    Kl , K2 , K3 , L , NR , NV , NDRS , NDRB , I Z , MASS , 
&    NRDOT,NVDOT 
DIMENSION    A(4,4),FVl(4),IVl(4),ZZZ(4,4) ,WR(4) ,WI(4) 

OPEN    (11, FILE='BIFl .RES' , STATUS- 'NEW ) 

Vehicle    Parameters: 
WEIGHT=435.0 


IZ 

=  45.0 

L 

=  7.3 

RHO 

=  1.94 

G 

=  32.2 

XG 

=0.0104 

MASS 

=WEIGHT/G 

YRDOT 

=-0.00178*0 

YVDOT 

=-0.03430*0 

YR 

=+0.01187*0 

YV 

=-0.03896*0 

YDRS 

=+0.02345*0 

YDRB 

=+0.02345*0 

NRDOT 

=-0.00047*0 

NVDOT 

=-0.00178*0 

NR 

=-0.01022*0 

NV 

=-0.00769*0 

NDRS 

=-0.337*0.0 

NDRB 

=+0.283*0.0 

5*RHO*L**4 

5*RHO*L**3 

5*RHO*L**3 

5*RHO*L**2 

5*RHO*L**2 

5*RHO*L**2 

5*RHO*L**5 

5*RHO*L**4 

5*RHO*L**4 

5*RHO*L**3 
02345*0. 5*RHO*L**3 
02  34  5*0 . 5*RHO*L**3 


WRITE  (*,1001) 

READ  (*,*)      TCMIN,TCMAX, ITC 

WRITE  (*,1003) 

READ  (*,*)     U 


DH   =( 

AA11=( 

AA12=( 

AA21=( 

AA22=( 

BB11=( 

BB12=( 

BB21=( 

BB22=( 

IZ-NRDOT) * (MASS-YVDOT)- 

MASS*XG-YRDOT)*(MASS*XG-NVDOT) 

( IZ-NRDOT) *YV-(MASS*XG-YRDOT) *NV)/DH 

( IZ-NRDOT) *(-MASS+YR)- 

MASS*XG-YRDOT)*(-MASS*XG+NR) )/DH 

(MASS-YVDOT) *NV-(MASS*XG-NVDOT) *YV)/DH 

( MASS-YVDOT) * ( -MASS*XG+NR )- 

MASS*XG-NVDOT) *(-MASS+YR) )/DH 

( IZ-NRDOT) *YDRS-(MASS*XG-YRDOT) *NDRS)/DH 

( IZ-NRDOT) *YDRB-(MASS*XG-YRDOT) *NDRB)/DH 

( MASS-YVDOT) *NDRS-( MASS*XG-NVDOT) *YDRS )/DH 

( MASS-YVDOT )*NDRB-( MASS *XG-NVDOT)* YDRB )/DH 


WRITE  (*,1004) 

READ   (*,*)      RATIO 

BBl=BBll+RATIO*BBl2 
BB2=BB21+RATIO*BB2  2 


WRITE  (*,1005) 
READ   (*,*) 

EPS=1 .D-10 
ITMAX-1000 


TL 
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DO  1  1=1, I TC 

WRITE  (*,2001)  I,ITC 

TOTCMIN+(  I-l)*(TCMAX-TCMIN)/(  ITC-1  ) 

TCD=TC*L/U 

POLEC-1 .0/TCD 

ADl=3.0*POLEC 

AD2-3 .0*POLEC**2 

AD3=POLEC**3 

Al=BBl*U*U 

Bl=BB2*U*U 

C1--AD1- ( AA11+AA2  2 ) *U 

A2=(BB1*AA2  2-BB2*AA12)*U**3 

B2-(BB2*AAll-BBl*AA21 )*U**3 

Kl=AD3/( (BB2*AAll-BBl*AA21 )*U**3) 

C2=AD2-( AA11*AA2  2-AA12*AA21 )*U**2+BB2*U*U*Kl 

K2=(Cl*B2-C2*Bl )/( Al *B2-A2 *Bl ) 

K3=(C2*A1-C1*A2)/(A1*B2-A2*B1 ) 
C 

Al-  (AAll*BB2-AA21*BBl )*U*U*U*Kl 

A2=  (AA11*AA2  2-AA12*AA21 ) *U*U+ ( AAl 1 *BB2-AA21 *BBl )*U*U*U*K3+ 
&        (AA2  2*BBl-AAl2*BB2)*U*U*U*K2-BB2*U*U*Kl 

A3=-(AAll+AA22)*U-(BBl*K2+BB2*K3)*U*U 

BO-  ( AAll*BB2-AA21*BBl ) *U*U*U*U*Kl 

Bl=-( AA12*BB2-AA22*BB1 )*U*U*U*Kl-BB2*U*U*U*Kl 

B2=-BB1*U*U*K1 
C 

C         COMPUTE  INITIAL  APPROXIMATION  FOR  OMEGA  (TL-0) 
C 

AQ=B2*A3-Bl 

BQ=B1*A2-B2*A1-B0*A3 

CQ-B0*Al 

DET-BQ**2-4 .D0*AQ*CQ 

IF  (DET. LT. 0.D0)  GO  TO  1 

ZTl=(-BQ+DSQRT(DET) )/(2.0D0*AQ) 

ZT2=(-BQ-DSQRT( DET) )/( 2 . 0D0*AQ) 

IF  (ZTl .GE.O.DO)  OM-DSQRT( ZT1 ) 

IF  (ZT2 .GE. 0.D0)  OM=DSQRT( ZT2 ) 

ALPHAl-AQ 

ALPHA2-BQ 

ALPHA3-CQ 

BETA1--B2 

BETA2=B0+B2*A2-B1*A3 

BETA3=Bl*Al-B0*A2 
C 

C        COMPUTE  EXACT  OMEGA  USING  NEWTON'S  ITERATION 
C 

OMOLD=OM 
3    F=(BETAl*OMOLD**5+BETA2*OMOLD**3+BETA3*OMOLD) *DSIN( OMOLD*TL) 
&    +(ALPHAl*OMOLD**4+ALPHA2*OMOLD**2+ALPHA3 ) * DCOS ( OMOLD*TL) 

FPRIME=( 5 . D0*BETAl*OMOLD**4+3 . DO  * BETA2 *OMOLD* * 2  +  BETA3 ) 
&    *DSIN(OMOLD*TL)  +  ( BETAl *OMOLD* * 5  +  BETA2 *OMOLD* *  3  + BETA 3 *OMOLD ) 
&    *TL*DCOS(OMOLD*TL)  +  ( 4 . D0*ALPHA1 *OMOLD* *  3  +  2 . DO *ALPHA2 *OMOLD ) 
&    *DCOS(OMOLD*TL)-( ALPHAl *OMOLD* * 4+ALPHA2 *OMOLD* * 2+ALPHA3 ) 
&    *TL*DSIN(OMOLD*TL) 

IF  ( FPRIME.EQ.O.DO)  STOP  1112 

IF  ( F. EQ.O .DO )  GO  TO  2 

OMNEW=OMOLD-F/FPRIME 

OMDI F-DABS ( OMNEW-OMOLD ) 

IF  (OMDIF.LE.EPS)  GO  TO  2 
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OMOLD-OMNEW 

IT-IT+1 

IF  ( IT.GT. ITMAX)  STOP  1111 

GO  TO  3 

2    OM-OMNEW 

XDNUM-DSQRT( (B0-B2*OM**2)**2+(B1*OM)**2) 
XDDEN=OM*DSQRT(  ( A1-A3 *OM* *  2 ) * *2+ ( A2 *OM-OM* *  3 ) *  *  2 ) 
XD=XDNUM/XDDEN 
XD=XD/L 

WRITEdl,*)  XD,TC,OM 

1  CONTINUE 

1001  FORMAT  ('  ENTER  MIN,  MAX,  AND  INCREMENTS  OF  TC ' ) 

1003  FORMAT  ('  ENTER  U') 

1004  FORMAT  ('  ENTER  BOW/STERN  RUDDER  RATIO') 

1005  FORMAT  ('  ENTER  TIME  LAG') 
2001  FORMAT  (215) 

END 
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